In [ ]:
import numpy as np
import pandas as pd
from tqdm import tqdm
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "notebook+pdf"
System Variables
In [ ]:
# Event codes
ARRIVAL=0
DEPARTURE=1
event_codes = {ARRIVAL: "ARRIVAL", DEPARTURE: "DEPARTURE"}
# Patient status codes
PENDING=0
ADMITTED=1
RELOCATED=2
REJECTED=3
DISCHARGED=4
status_codes = {PENDING: "PENDING", ADMITTED: "ADMITTED", RELOCATED: "RELOCATED", REJECTED: "REJECTED"}
# Simulation parameters
SIM_TIME = 31 # Simulation time in days
# Ward parameters
WARDS = dict(
names = ['A', 'B', 'C', 'D', 'E', 'F'],
lams = [14.5, 11.0, 8.0, 6.5, 5.0, 13.0],
mu_invs = [2.9, 4.0, 4.5, 1.4, 3.9, 2.2],
urgency_points = [7, 5, 2, 10, 5, 0],
)
NUM_WARDS = len(WARDS['names'])
RELOCATION_PROBABILITIES = np.array([
[0.0, 0.05, 0.10, 0.05, 0.80, 0.00],
[0.2, 0, 0.50, 0.15, 0.15, 0.00],
[0.30, 0.20, 0, 0.20, 0.30, 0.00],
[0.35, 0.30, 0.05, 0, 0.3, 0.00],
[0.20, 0.10, 0.60 ,0.10, 0, 0.00],
[0.20, 0.20, 0.20, 0.20, 0.20 ,0]
])
Penalty Computation
In [ ]:
def penalty(urgency_points, counts):
rates = (counts[RELOCATED] + counts[REJECTED]) / (counts[ADMITTED] + counts[RELOCATED] + counts[REJECTED])
return np.dot(rates, urgency_points)
Simulation Setup
In [ ]:
def simulation_setup(
sim_time=SIM_TIME,
arrival_sampler=lambda lam: np.random.exponential(1/lam),
los_sampler=lambda mu_inv: np.random.exponential(mu_inv)
):
# Patients
patients = dict( ID=[], type=[], ward=[], status=[], burn_in=[] )
# Precompute arrival and departure times for each ward
events = dict( event_type=[], patient_ID=[], time=[] )
patient_id = 0
for i, (lam, mu_inv) in enumerate(zip(WARDS['lams'], WARDS['mu_invs'])):
# Sample patients for ward
clock = 0
while clock < sim_time:
# Sample arrival time
clock += arrival_sampler(lam)
# Add patient to event list if arrival time is before simulation end
if clock <= sim_time:
# Event
events['time'] += [clock, clock+los_sampler(mu_inv)]
events['patient_ID'] += [patient_id]*2
events['event_type'] += [ARRIVAL, DEPARTURE]
# Patient
patients['ID'].append(patient_id)
patients['type'].append(i)
patients['ward'].append(i)
patients['status'].append(PENDING)
patients['burn_in'].append(0)
# Update patient ID
patient_id += 1
# Sort events by time
idx = np.argsort(events['time'])
for key in events.keys():
events[key] = [events[key][i] for i in idx]
return patients, events
Inspect the Event List
In [ ]:
# Setup simulation
patients, events = simulation_setup()
# Create dataframes and map codes to names
patient_df = pd.DataFrame(patients); patient_df['status'] = patient_df['status'].map(status_codes)
patient_df['type'] = patient_df['type'].map({i: WARDS['names'][i] for i in range(NUM_WARDS)})
event_df = pd.DataFrame(events); event_df['event_type'] = event_df['event_type'].map(event_codes)
# Print dataframes
print("Patients:")
print(patient_df.head().to_string(index=False), "\n")
print("Events:")
print(event_df.head().to_string(), "\n")
Patients: ID type ward status burn_in 0 A 0 PENDING 0 1 A 0 PENDING 0 2 A 0 PENDING 0 3 A 0 PENDING 0 4 A 0 PENDING 0 Events: event_type patient_ID time 0 ARRIVAL 0 0.011755 1 ARRIVAL 432 0.017686 2 ARRIVAL 1074 0.045629 3 DEPARTURE 1074 0.049026 4 ARRIVAL 794 0.069206
Simulation Loop
In [ ]:
def simulate(patients, events, capacities, burn_in_period=15):
# Clear patient status
N = len(patients['ID'])
patients['ward'] = [None]*N
patients['status'] = [PENDING]*N
# List for saving ward states
states = []
# Dict for saving counts
counts = {
ADMITTED: np.zeros(NUM_WARDS, dtype=int),
RELOCATED: np.zeros(NUM_WARDS, dtype=int),
REJECTED: np.zeros(NUM_WARDS, dtype=int),
}
# Simulate loop
occupancy = np.zeros(NUM_WARDS, dtype=int)
for i, (event_type, patient_id) in enumerate(zip(events['event_type'], events['patient_ID'])):
# Check if patient is in burn-in period
if events['time'][i] < burn_in_period:
patients['burn_in'][patient_id] = True
# Get patient type
patient_type = patients['type'][patient_id]
if event_type == ARRIVAL:
# Admit patient
if occupancy[patient_type] < capacities[patient_type]:
ward = patients['type'][patient_id]
status = ADMITTED
else:
# Relocate patient (or reject if alternative ward is full)
ward = np.random.choice(NUM_WARDS, p=RELOCATION_PROBABILITIES[patient_type])
status = RELOCATED if occupancy[ward] < capacities[ward] else REJECTED
# Update patient
if status != REJECTED:
occupancy[ward] += 1
patients['ward'][patient_id] = ward
patients['status'][patient_id] = status
if events['time'][i] >= burn_in_period:
counts[status][patient_type] += 1
# Discharge patient
elif patients['status'][patient_id] != REJECTED:
occupancy[patients['ward'][patient_id]] -= 1
states.append(occupancy.copy())
events['states'] = states
# Check if simulation ended with non-zero occupancy
if any(occupancy != 0):
print('Error: Simulation ended with non-zero occupancy')
return dict(states=states, counts=counts, penalty=penalty(WARDS['urgency_points'], counts))
Run Simulation
In [ ]:
patients, events = simulation_setup()
capacities = [49, 28, 22, 16, 16, 34]
sim_out = simulate(patients, events, capacities)
states = sim_out['states']
fig = go.Figure()
for i, state in enumerate(np.array(events['states']).T):
fig.add_trace(go.Scatter(x=events['time'], y=state, mode='lines', name=WARDS['names'][i]))
fig.update_layout(title='Ward occupancies', xaxis_title='Time', yaxis_title='Occupancy', width=800)
fig.show()
Gradient Computation
In [ ]:
def compute_gradients(capacities):
M = len(capacities)
grads = np.zeros(M)
n = 10
def simulate_penalty(capacity_adjustment):
total_penalty = 0
for _ in range(n):
capacities[i] += capacity_adjustment
patients, events = simulation_setup()
sim_out = simulate(patients, events, capacities)
total_penalty += sim_out['penalty']
capacities[i] -= capacity_adjustment
return total_penalty / n
for i in range(M):
p_small = simulate_penalty(-1)
p_large = simulate_penalty(1)
grads[i] = (p_large - p_small) / 2
return grads
def gradient_descent(capacities, num_beds_F=34):
beds_to_move = num_beds_F - capacities[-1]
penalties = np.zeros(beds_to_move+1)
F_rel_prob = np.zeros(beds_to_move+1)
# Use the samee simulations to estimate performance metrics
sim_setup_samples = [simulation_setup() for _ in range(100)]
# Realocate beds to F
for i in tqdm(range(beds_to_move+1)):
# Simulate and compute performance metrics
sim_out = [simulate(patients, events, capacities) for patients, events in sim_setup_samples]
penalties[i] = np.mean(np.array([out['penalty'] for out in sim_out]))
num = [out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
den = [out['counts'][ADMITTED][-1] + out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
F_rel_prob[i] = np.mean([num/den for num, den in zip(num, den)])
# Compute gradients and update capacities
grads = compute_gradients(capacities)
idx = np.argmax(grads[:-1])
capacities[idx] -= 1
capacities[-1] += 1
# Correct for last iteration
capacities[idx] += 1; capacities[-1] -= 1
return dict(capacities=capacities, penalties=penalties, F_rel_prob=F_rel_prob)
Optimizing Bed Allocation
In [ ]:
# Run gradient descent for different bed totals (Takes a few minutes)
# out165 = gradient_descent([55, 40, 30, 20, 20, 0])
# capacities165 = out165['capacities']
# out170 = gradient_descent([55, 40, 30, 20, 20, 5])
# capacities170 = out170['capacities']
# out180 = gradient_descent([55, 40, 30, 20, 20, 15])
# capacities180 = out180['capacities']
In [ ]:
out165 = {'capacities': [47, 33, 22, 15, 14, 34],
'penalties': np.array([7.49802066, 7.55358105, 7.58117201, 7.59598008, 7.59898524,
7.62337703, 7.66993802, 7.67370977, 7.81047128, 7.89989524,
7.9138084 , 7.8682587 , 8.11521142, 8.06129237, 8.05570368,
8.14991408, 8.34454141, 8.44293029, 8.44261769, 8.51071874,
8.58172614, 8.50190952, 8.68852844, 8.85626831, 8.8376133 ,
8.75429536, 9.10895572, 9.24257824, 9.1926232 , 9.26930196,
9.35233008, 9.45271561, 9.37556658, 9.53496749, 9.97145061]),
'F_rel_prob': np.array([1. , 0.96736859, 0.93616246, 0.90050175, 0.86772832,
0.83351419, 0.80002296, 0.7673059 , 0.73387623, 0.70163417,
0.66996514, 0.63766473, 0.60441516, 0.57166332, 0.54084904,
0.50704336, 0.47572081, 0.44561348, 0.41672475, 0.38703384,
0.35834503, 0.32898252, 0.30078676, 0.27299778, 0.24596464,
0.21982852, 0.19390895, 0.17003065, 0.14740803, 0.12769683,
0.10878097, 0.09190907, 0.07592539, 0.06171736, 0.04995138])}
out170 = {'capacities': [48, 34, 22, 17, 15, 34],
'penalties': np.array([6.61726877, 6.60865111, 6.56312237, 6.69416874, 6.63901451,
6.77038572, 6.8701144 , 6.86576323, 6.91945093, 6.97832302,
6.93358418, 6.95258001, 6.91292444, 6.99449613, 6.93944662,
7.11645418, 7.11990933, 7.31588464, 7.20711269, 7.60083692,
7.58557202, 7.73770086, 7.72940449, 7.7162035 , 7.96621445,
8.06607762, 8.11213987, 8.29464205, 8.43521303, 8.51591592]),
'F_rel_prob': np.array([0.82980978, 0.7935083 , 0.76158163, 0.72749367, 0.69373752,
0.66053872, 0.62795306, 0.59538922, 0.56263003, 0.53105079,
0.50077906, 0.4688307 , 0.43661128, 0.40846731, 0.38075466,
0.34940134, 0.31830065, 0.28883195, 0.26118983, 0.23610832,
0.21077812, 0.18653657, 0.16454685, 0.14484725, 0.12425673,
0.10599562, 0.08915169, 0.07332146, 0.05900622, 0.04685821])}
out180 = {'capacities': [52, 35, 25, 17, 17, 34],
'penalties': np.array([5.59263092, 5.57043543, 5.61606463, 5.7193527 , 5.67450016,
5.80707938, 5.70256907, 5.78518861, 5.93776501, 6.06845177,
6.14692932, 6.13932206, 6.30083232, 6.33349891, 6.37718681,
6.612897 , 6.76277317, 6.89640734, 7.18052889, 7.12720848]),
'F_rel_prob': np.array([0.50123168, 0.46873188, 0.43965553, 0.40936768, 0.37708796,
0.34639364, 0.31709197, 0.28659546, 0.25933879, 0.23323819,
0.2071624 , 0.18398445, 0.16063215, 0.1391825 , 0.11869791,
0.10052395, 0.08263447, 0.06724811, 0.05403682, 0.04260908])}
Theoretical Optimal Number of Beds in F
In [ ]:
def erlangB_formula(m):
if m == 0:
return 1
lam = 13; s = 2.2
A = lam*s
def power_div_factorial(i):
return np.exp(i*np.log(A) - np.sum(np.log(range(1, i+1))))
return power_div_factorial(m)/(np.sum([power_div_factorial(i) for i in range(m+1)]))
In [ ]:
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=[erlangB_formula(m) for m in range(35)], mode='lines', name='Erlang B'))
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['F_rel_prob'], mode='lines', name='Simulation'))
fig.update_layout(title='Probability of relocating in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Probability', width=800)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['penalties'], mode='lines', name='Simulation'))
fig.update_layout(title='Penalty as a function of number of beds in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Penalty', width=800)
fig.show()
Estimate Rates for Admission, Relocation and Rejection
In [ ]:
def estimate_rates(samples, capacities):
num_samples = len(samples)
rates = np.zeros((num_samples, 3, 6))
for i in range(num_samples):
patients, events = samples[i]
sim_out = simulate(patients, events, capacities)
counts = sim_out['counts']
rates[i, 0] = counts[ADMITTED]
rates[i, 1] = counts[RELOCATED]
rates[i, 2] = counts[REJECTED]
rates /= rates.sum(axis=(1, 2))[:, None, None]
rates = rates.mean(axis=0)
rates = pd.DataFrame(rates.T, index=WARDS['names'], columns=['ADMITTED', 'RELOCATED', 'REJECTED'])
return rates
In [ ]:
num_samples = 100
sim_setup_samples = [simulation_setup() for _ in range(num_samples)]
for n, cap in zip([165, 170, 180], [out165['capacities'], out170['capacities'], out180['capacities']]):
rates = estimate_rates(sim_setup_samples, cap)
print(f'{n} beds:')
print(rates.to_string(), '\n')
165 beds: ADMITTED RELOCATED REJECTED A 0.198867 0.016467 0.034129 B 0.118000 0.035187 0.035894 C 0.058521 0.047748 0.031763 D 0.084323 0.016815 0.012998 E 0.035033 0.025906 0.024932 F 0.213094 0.005663 0.004661 170 beds: ADMITTED RELOCATED REJECTED A 0.203240 0.015898 0.030325 B 0.122021 0.034313 0.032746 C 0.058817 0.050653 0.028562 D 0.093151 0.012455 0.008530 E 0.037995 0.025517 0.022359 F 0.213094 0.006083 0.004241 180 beds: ADMITTED RELOCATED REJECTED A 0.220776 0.011276 0.017411 B 0.127825 0.035854 0.025402 C 0.069993 0.047934 0.020105 D 0.097713 0.010247 0.006176 E 0.046178 0.023315 0.016378 F 0.213094 0.006721 0.003604
Variance and Distribution of LOS
In [ ]:
var_scales = np.linspace(1, 10, 10)
penalties = np.zeros(10)
mu_invs = WARDS['mu_invs']
for i, var_scale in enumerate(var_scales):
WARDS['mu_invs'] = [[np.log(mu_inv), np.sqrt(var_scale/mu_inv**2)] for mu_inv in mu_invs]
los_sampler = lambda mu_inv: np.random.lognormal(mean = mu_inv[0], sigma = mu_inv[1])
sim_setup_samples = [simulation_setup(los_sampler=los_sampler) for _ in range(100)]
penalties[i] = np.mean([simulate(patients, events, capacities)['penalty'] for patients, events in sim_setup_samples])
WARDS['mu_invs'] = mu_invs
fig = go.Figure()
fig.add_trace(go.Scatter(x=var_scales, y=penalties, mode='lines'))
fig.update_layout(title='Penalty as a function of variance scale', xaxis_title='Variance scale', yaxis_title='Penalty', width=800)
fig.show()